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Abstract 

We study numerically the spatial properties of relaxation clusters in a two dimensional 
sandpile automaton with dynamic rules depending stochastically on a parameter p, which 
models the effects of static friction. In the limiting cases p = 1 and p = the model 
reduces to the critical height model and critical slope model, respectively. At p = p c , a 
continuous phase transition occurs to the state characterized by a nonzero average slope. 
Our analysis reveals that the loss of finite average slope at the transition is accompanied 
by the loss of fractal properties of the relaxation clusters. 
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1 Introduction 

The term self-organized criticality (SOC) (Bak 1987, Bak 1988) applies to certain nonlin- 
ear extended dynamical systems which tend to self-organize into a steady state with long 
range spatial correlations, analogous to the critical state close to equilibrium second-order 
phase transitions. The dynamics of these systems consists of series of events (avalanches) 
in which the system is repeatedly perturbed and let to relax according to certain micro- 
scopic dynamic rules. After a certain large number of time steps, the system " learns" its 
response to the external perturbation. Clusters of relaxed sites have been shown to have 
fractal properties (Janosi 1994). Therefore, the study of spatial properties of relaxation 
clusters offers an additional possibility to analyze the character of the dynamical process 
itself. 

Sandpile automata models (Bak 1987, Dhar 1989, Manna 1991) are well known pro- 
totype models exhibiting SOC. By adding a particle from the outside, a sandpile is 
perturbed and the perturbation may lead to instabilities at neighboring sites if the local 
height of sand exceeds a critical value h c — critical height model (CHM) (Bak 1987, 
Dhar 1989), or if the local slopes exceed a critical value a c — critical slope model (CSM) 
(Liibeck 1993, Manna 1991). The self-tuning of the system in the case of the CHM is 
accompanied with a state of a zero average slope. A finite average slope characterizes a 
CSM. 

In the present work we study a two dimensional sandpile-type automaton with pre- 
ferred direction in which the updating rules are tuned continuously by changing a pa- 
rameter p between CHM (for p = 1) and CSM (for p — 0). For all values < p < 1 
the dynamics is stochastic. According to these dynamic rules, the relaxation clusters 
may have holes and dendritic forms. We present the results of numerical investigations 
of the spatial properties, i.e. distributions of size and length, and fractal dimension of 
the relaxation clusters. Our conclusion is that there is a treshold value of the probability 
parameter, p = p c , above which the system behaves like a CHM, characterized by zero 
net slope. However, for p < p c the dynamics leads to the state with a nonzero net slope. 
The phase transition between these two states of the automaton is characterized by a 
continuous appearance of the average slope at p < p c . At the same time, on approaching 
p = p c from above, the fractality of the relaxation clusters disappears. 
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2 Model and Numerical Simulations 



We consider a two-dimensional sandpile model on a square lattice of size L x L and 
integer variables h(i,j), representing the local height. We assume a directed dynamics, 
i.e. particles are restricted to flow in the downward direction (increasing i). According 
to the widespread 'sandpile language' the first row (i = 1) and the last row (i = L) 
represent the top and the bottom of the pile, respectively. To minimize the influence of 
the horizontal boundaries we limit our investigations to periodic boundary conditions in 
this horizontal direction (j-direction). Any site of the lattice has two downward and two 
upward next neighbours, namely 

h(i + l,j±) and h(i-l,j ± ) with J± =j± {l±{ -~ l)t \ (1) 

We perturb the system by adding particles at a random place on the top of the pile 
according to 

h{\,j) i— > h(\,j) + 1 , with random j. (2) 

A site is called unstable if the height h(i,j) or at least one of the two slopes 

<r(i, 3±) = Kh j) ~ h(i + 1, j±) (3) 

exceeds a critical value, i.e if h(i,j) > h c or <r(i,j±) > cr c , respectively. 

In the case of the critical slope condition toppling takes place until both slopes become 
subcritical. In this toppling process particles drop alternative to the downward next 
neighbours if both slopes are unstable. If only one slope exceeds a c particles drop to the 
corresponding downward neighbour. 

In contrast to the critical slope the critical height conditions has a stochastic character. 
If the local height exceeds the critical value h c toppling occurs only with the probability 
p, and then two particles drop to the two downward neighbours. 

Because each relaxing cell changes the heights and slopes of its four next neighbouring 
sites the stability conditions are applied at these four 'activated' sites in the next updating 
step. Toppling may take place after adding a particle on the first row. We first apply 
the slope and then the height stability condition in parallel for each activated site. An 
avalanche stops if all activated sites are stable. Then we start again according to Eq. (2). 
One can interpret 1 — p as being due to static friction between the sand grains, which 
prevents toppling even if the height exceeds the critical value. In the case p — 1, our 
model is identical to the model of Dhar and Ramaswamy (Dhar 1989), which exhibits a 
robust SOC behavior. In the limit p — 0, our directed rules lead to a steady state where 
all slopes are equal, i.e., <r(i,j±) — a c — 1. Thus, if a particle is added at the first row it 
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performs a directed random walk downward the pile until it reaches the boundary and 
drops out of the system — no SOC can occur in this limit. 

We fix the critical height h c and the critical slope a c and restrict ourselves to h c = 2 
and a c = 8, although it should be emphasized that the results depend on the choice of 
these parameters. 




Figure 1 Snapshots of sandpiles for p — 1 (upper left), 0.7 (upper right), 0.3 (lower 
left), and 0.1. In all cases the upper left boundary corresponds to the first row of the pile 
(i = 1). 

Starting from an empty lattice of linear size L we add approximately L 3 particles to 
equilibrate the system before taking measurements. In each case the stationarity of the 
system was checked by determining that the average height <h(i,j)> has reached a con- 
stant value. 
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3 Results 



In Figure 1 we show four snapshots of the sandpile for different values of p. The upper 
left picture represents the pile in a steady state of the pure CHM (p = 1). Only the 
heights and 1 are present in the critical state. With decreasing p < 1 a certain number 
of heights with values h(i,j) > h c remain in the interior of an avalanche leading to a 
rough surface of the pile (see e.g. the case p = 0.7 of Figure 1). 

For still lower values of p the roughness becomes more dramatic and eventually for 
p = 0.3 (lower left) the surface shows large fluctuations, indicating a change of the system 
behavior. In the case p — 0.1, displayed in the lower right picture, the heights grow up to 
the full size required by the CSM rules with fluctuations around the well defined average 
slope. We now analyze this behavior quantitatively by considering the average height 
(h(i)) as a function of distance from the top row i, defined as 

(h(i)) = L-^WJ)), (4) 

where (h(i,j)) is determined as the average over total number of time steps. In the 
interior of the pile the average height is independent of % for all values of p > 0.3, i.e. 
(h(i)) = (h), where (h) depends on the probability parameter p. The normalized average 
heights are plotted in Figure 2a vs. the distance %. Except of the close vicinity of the 
boundaries the curves corresponding to different values of p > 0.3 collapse into a single 
curve (h(i)) / (h) = 1. 

However, this is no longer the case for p < 0.3, where the deviations from (h) char- 
acteristic for the behavior at the boundaries proliferates into the interior of the pile. 
Instead of constant values of the average heights we find for p < 0.3 constant values of 
the average slopes, which are defined as a(i) = (h(i)) — (h(i + 1)). In Figure 2b the 
normalized average slopes are shown. Except of the deviations due to the boundaries all 
curves more or less collapse to a(i)/a(i = 50) = 1. Simulations for some values of p and 
different automata sizes (L = 50, 100, 200) show that these deviations are independent of 
L, i.e. they are true boundary effects. In the interior part of the automaton the average 
slopes are independent of %. 

For p < 0.3, the slopes obey a translation invariance similar to that of the heights in 
the opposite limit. The average slopes become a function of i if the probability parameter 
p is close to p — 0.3. 
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Figure 2 (a) Normalized average height (h(i)) at distance % from the top of the pile 
plotted vs. the distance % for various values of p between 0.9 (flat curve) and 0.295 (curve 
with profile), (b) Normalized average slope (cr(i)) for various values of p between 0.1 (flat 
curve) and 0.29 (curve with profile). 



Since a zero net slope characterizes a CHM and a non zero net slope a CSM, respec- 
tively, we interpret the average slopes 
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i=20 



as an order parameter, which describes the transition from the critical height to the 
critical slope regime. In this definition we cut off the first and last 20 values in order to 
minimize the influence of the boundaries. The results for (a) which are obtained in this 
way are independent of L as additional simulations for L = 200 have shown. Of course 
it is also irrelevant if the cut off length is larger than 20. 
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The true nature of the transition as well as the transition point p c should be determined 
from the p-dependence of the average slope itself in the vicinity of p c . We determined 
the average net slope according to Eq. (||) for different values of p. 

The results are shown in Figure 3. As the average net slope plays the role of an order 
parameter in our model, we see from Figure 3 that a continuous transition occurs between 
the two regimes of the system. Assuming the following form 

(v)~(Pc-pT, (6) 

we determine the transition point p c = 0.293 ± 0.002 and the exponent a = 0.8 ± 0.05 
from the data in Figure 3. 




Figure 3 Average slope of the pile (a) plotted vs. probability p. The solid line of the 
inset corresponds to a fit according to equation with p c = 0.293 and a = 0.8. 

Next we consider the spatial structure of the avalanches. One relaxation cluster con- 
sists of all sites that topple during one event following one particle being added at the top 
row. In Figure 4 we show six representative snapshots of the relaxation clusters taken 
from simulations for different values of p. For p = 1 the avalanches are compact. With 
decreasing p some supercritical sites may remain in the interior of the avalanche due 
to the stochastic character of the dynamic rules, leading to holes and branching of the 
clusters. For values of p which are lower than p c « 0.3 this structure of the avalanches is 
lost. In this region the net nonzero slope is developed and the probability of triggering 
an after-avalanche increases dramatically. The system relaxes occasionally through huge 
avalanches involving almost all sites. In the limiting case of p = the system shows 
rather simple behavior, since the rules of the CSM are deterministic and fluctuations 
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arrond the average slope a c — l are absent. Each added particle performs a random walk 
down the whole automaton. 
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Figure 4 Snapshots of relaxation clusters for a few values of p for a linear system size 
of L = 80. Sites which toppled during an avalanche are marked as black. The base of 
every picture represents the bottom of the pile (i = L). 

We next analyze the scaling properties of the clusters, which we expect to reflect these 
observations. We calculate the average size of all clusters of length /, (s)(l), which scales 
as a power of the length /, i.e., 

(s)(l) ~ h . (7) 

We emphasize that this definition of d\\ describes only the scaling along the prefered 
direction. Actually the avalanches in our model are self-affine fractals, i.e. the fractal 
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dimension along the perpendicular direction d± is not equal to d\\ . A detailed discussion 
of this anisotropic behavior will be published elsewhere. 
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Figure 5 (a) Double logarithmic plot of the average size of clusters (s)(l) vs. length 
/ for various values of p as indicated, (b) Parallel fractal dimension du of the relaxation 
clusters plotted vs. probability p. For p > 0.33 the error-bars are smaller than the sym- 
bols. 

The exponent in Eq. fl?]), which determines the scaling properties of relaxation clus- 
ters, depends on the parameter p (see Figure 5). We examined this p-dependence of the 
fractal dimension intensively, the obtained results are shown in Figure 5. In agreement 
with the pictures of the avalanche shapes in Figure 4 du increases with decreasing p. 
From d|| = 1.469 for p = 1, which is in a good agreement with the known value d\\ = | 
for the pure directed CHM (Tadic 1992), the fractal dimension grows with decreasing p 
and reaches its maximum value dn = 2 for p < 0.3. In this way the critical height regime 



9 



is accompanied by a self-afnne fractal structure of the relaxation cluster. This fractal 
behavior vanishes below the transition point p c where the critical slope condition rules 
the dynamics. 

4 Conclusions 

In conclusion, we have shown that our stochastic sandpile automaton exhibits a continu- 
ous phase transition in its steady state, as the probability parameter p is varied through 
the transition point p c = 0.293. The steady states for p < p c are characterized by a finite 
net average slope. For p > p c the average slope remains zero. Fractality of the relaxation 
clusters, measured by the deviation of their fractal dimension d\\ from the Euclidean di- 
mension d — 2, is shown to vanish at the transition point, strongly suggesting that the 
fractal character of the relaxation clusters is closely related to the CHM regime. 

The analysis of the avalanche distributions and especially the question whether the 
model displays SOC behavior for p < 1 remains beyond the scope of the present paper 
and will be published elsewhere. 

Acknowledgments 

This work was supported in part by the Deutsche Forschungsgemeinschaft through Son- 
derforschungsbereich 166, Germany, and by the Ministry of Science and Technology of 
the Republic of Slovenia. We would like to thank U. Nowak for critical reading the 
manuscript. One of us (B.T.) wishes to acknowledge the hospitality at the University of 
Duisburg, where this work was initiated. 



10 



References 

Bak P., Tang C. and Wiesenfeld K. (1987) Self- Organized Criticality: An Explanation of 

1/f Noise, Physical Review Letters, 59, 381. 
Bak P., Tang C. and Wiesenfeld K. (1988) Self-organized criticality, 

Physical Review A, 38, 364. 
Dhar D. and Ramaswamy R. (1989) Exactly Solved Model of Self-Organized Critical 

Phenomena, Physical Review Letters, 63, 1659. 
Janosi I.M. and Csirok A. (1994) Fractal clusters and self-organized criticality, 
Fractals, 2, 153. 

Liibeck S., Usadel K.D. (1993) SOC in a class of sandpile models with stochastic 

dynamics, Fractals 1, 1030. 
Manna S. S. (1991) Critical exponents of the sand pile models in two dimensions, 

Physica A 179, 249. 
Tadic B., Nowak U., Usadel K.D., Ramaswami R., and Padlewski S. (1992) 

Scaling behavior in disordered sandpile automata, Physical Review A, 45, 8536. 



11 



